//
//  Simulation.cpp
//  phy_four
//
//  Created by Snir Gazit on 3/12/12.
//  Copyright (c) 2012 __MyCompanyName__. All rights reserved.
//

#include "Simulation.h"

typedef boost::chrono::system_clock::time_point TimePoint;
typedef boost::chrono::duration<double, boost::ratio<3600>> Hour;
typedef boost::chrono::duration<double, boost::micro> MicroSec;

double calcWR(double mu, double U, int m) {
	double c1, c2, c3, md = double(m), arg;
	arg = 0.25 * pow(mu, 2) / U;
	c1 = 0.25 * pow(U, -(5. / 4.) - double(md) / 2);
	c2 = -mu * gsl_sf_gamma((5. / 4.) + md / 2)
			* gsl_sf_hyperg_1F1((5. / 4.) + md / 2, 3. / 2., arg);
	c3 = sqrt(U) * gsl_sf_gamma((3. / 4.) + md / 2)
			* gsl_sf_hyperg_1F1((3. / 4.) + md / 2, 1. / 2., arg);
	return c1 * (c2 + c3);
}

void calcWeight(double muIn, double UIn, double g, int maxBond,
		ostream &outFile, BoltzWeightType &Rweight) {
	double mu = muIn * g, U = UIn * g;
	double n, c, w;
	c = calcWR(mu, U, 0);
	for (int i = 1; i < MAX_BOND; i++) {
		n = calcWR(mu, U, i);
		w = (n / c);
		c = n;
		Rweight[i - 1] = w;
	}
}
// Worm Alg

void OnModel::init(QmcParameters &qmcParameters) {
	//Init Params
	qmcParameters_ = qmcParameters;
	string runNum = qmcParameters_.runNum.getVal();
	// OutFile
	string fileName = qmcParameters_.outPath.getVal() + "/out" + runNum
			+ ".txt";
	outFile_.open(fileName.c_str());
	outFile_ << "Starting N=3 Simulation, Run Number - " + runNum << endl;
	// Calc Weights
	qmcParameters_.g.setVal(1. / (qmcParameters_.g.getVal()));
	if (qmcParameters_.simType.getVal() == "HEISENBERG")
		boltzWeight_.fill(1.);
	else
		calcWeight(qmcParameters_.mu.getVal(), qmcParameters_.U.getVal(),
				qmcParameters_.g.getVal(), qmcParameters_.maximalBonds.getVal(),
				outFile_, boltzWeight_);
	boost::array<double, 2> prob = { 2./3.,1./3.};
	// Init Lattice
	boost::array<int, DIM> sizes = qmcParameters_.Llist.getVal();
	lat_.init(sizes);
	lat_.setSectorType(ZSECTOR);
	// Init Rng
	eng_.seed(qmcParameters_.seed.getVal());
	// Init Moves
	//	createMoveXY_.init(&eng_, &lat_, qmcParameters_, XYWeight);
	//	createMoveIsing_.init(&eng_, &lat_, qmcParameters_, IsingWeight);
	//  removeMoveXY_.init(&eng_, &lat_, qmcParameters_, XYWeight);
	//removeMoveIsing_.init(&eng_, &lat_, qmcParameters_, IsingWeight);
	createPairMoveIsing_.init(&eng_, &lat_, qmcParameters_, boltzWeight_);
	createPairMoveXY_.init(&eng_, &lat_, qmcParameters_, boltzWeight_);
	shiftMoveIsing_.init(&eng_, &lat_, qmcParameters_, boltzWeight_);
	shiftMoveXY_.init(&eng_, &lat_, qmcParameters_, boltzWeight_);
	chooseLoopType_.init(&eng_, DiscreteVariableParams(prob));
	// Initialize coupeling
	setCoupeling(qmcParameters_.g.getVal());
	// Measurements
	measures_.init(&lat_, qmcParameters_, &boltzWeight_, &outFile_);
	setCoupeling(qmcParameters_.g.getVal());
	// Read From File
	if (qmcParameters_.initFromFile.getVal())
		this->read();
	int totalNumberOfMeasure = (int) pow(
			(double) qmcParameters.binSize.getVal(),
			qmcParameters.numOfBins.getVal());
	//measures_.updateMinNumOfMeasure();
}
OnModel::~OnModel() {
	outFile_.close();
}

void OnModel::thermalize() {
	int numOfThermal = qmcParameters_.numOfThermal.getVal();
	acc CProb, XYProb;
	for (int i = 0; i < numOfThermal; i++) {
		nextMove();
		if (lat_.isClosed())
			CProb(1.);
		else
			CProb(0.);
		if (lat_.getSectorType() == GXY)
			XYProb(1.);
		else
			XYProb(0.);
	}
	bool conv;
	double cCProb = ba::mean(CProb), nCProb = 0.0;
	do {
		CProb = acc();
		XYProb = acc();
		for (int i = 0; i < numOfThermal; i++) {
			nextMove();
			if (lat_.isClosed())
				CProb(1.);
			else
				CProb(0.);
			if (lat_.getSectorType() == GXY)
				XYProb(1.);
			else
				XYProb(0.);
		}
		double nCProb = ba::mean(CProb);
		outFile_ << "Closed Prob " << nCProb << endl;
		outFile_ << "XY Prob " << ba::mean(XYProb) << endl;
		double ratio = fabs(1. - (cCProb / nCProb));
		conv = ratio < 0.01;
		cCProb = nCProb;
	} while (!(conv));
	//CreateLoopMove_.setAvWormLength(2 * ba::mean(wormLength));
	outFile_ << "Final Closed Prob - " << cCProb << endl;
	int closedRate = int(
			double(qmcParameters_.rateFactor.getVal())
					* double(lat_.getNumOfSites()) * cCProb);
	//	closedRate = 1;
	measures_.updateClosedRate(closedRate);
	outFile_ << "Closed Rate - " << closedRate << endl;
}

bool OnModel::run() {
	bool endTime, dumpFlag;
	namespace ba = boost::accumulators;
	ba::accumulator_set<double, ba::stats<ba::tag::mean> > wormLength;
	while (measures_.checkNotDone()) {
		nextMove();
		measures_.DoMeasurement(endTime, dumpFlag);
		if (dumpFlag) {
			this->dump(false);
		}
		if (endTime) {
			//this->dump(true);
			break;
		}
	}
	return !endTime;
}
void OnModel::nextMove() {
	SectorType cSectorType = lat_.getSectorType();
	if (cSectorType == GISING) {
		wormLength_ = 0;
		do {
			if (shiftMoveIsing_.makeMove())
				wormLength_++;
		} while (lat_.getSectorType() != ZSECTOR);
		wormLengthIsing_(wormLength_);
	} else if (cSectorType == GXY) {
		wormLength_ = 0;
		do {
			if (shiftMoveXY_.makeMove())
				wormLength_++;
		} while (lat_.getSectorType() != ZSECTOR);
		wormLengthXY_(wormLength_);
	} else {
		int sectorType = chooseLoopType_.draw();
		if (sectorType == 0) {
			if (createPairMoveXY_.makeMove())
				wormLength_++;
		} else {
			if (createPairMoveIsing_.makeMove())
				wormLength_++;
		}
	}
}
void OnModel::setCoupeling(double g) {
	calcWeight(qmcParameters_.mu.getVal(), qmcParameters_.U.getVal(), g,
			qmcParameters_.maximalBonds.getVal(), outFile_, boltzWeight_);
	shiftMoveXY_.setCoupeling(boltzWeight_, g);
	shiftMoveIsing_.setCoupeling(boltzWeight_, g);
	createPairMoveIsing_.setCoupeling(boltzWeight_, g);
	createPairMoveXY_.setCoupeling(boltzWeight_, g);
}

void OnModel::simulatedAnnealing() {
	outFile_ << " Startting Simulated Anneling " << endl;
	double gHighT = 1. / qmcParameters_.gHighT.getVal();
	int numOfSMA = qmcParameters_.numOfGSA.getVal();
	double dg = (qmcParameters_.g.getVal() - gHighT) / double(numOfSMA - 1);
	double nLengthIsing, nLengthXY, nTotal;
	int numOfThermal = qmcParameters_.numOfThermal.getVal();
	vector<double> gSA;
	int numIter = 0;
	for (int i = 0; i < qmcParameters_.numOfGSA.getVal(); i++)
		gSA.push_back(gHighT + dg * double(i));
	for (auto it_g = gSA.begin(); it_g != gSA.end();it_g++) {
		bool conv;
		double cg = (*it_g);
		setCoupeling(cg);
		outFile_ << "Current g -" << cg << endl;
		wormLengthXY_ = acc();
		wormLengthIsing_ = acc();
		int i = 0;
		while (i < numOfThermal) {
			nextMove();
			if (lat_.sectorType_ != ZSECTOR)
				i++;
		}
		nLengthXY = ba::mean(wormLengthXY_);
		nLengthIsing = ba::mean(wormLengthIsing_);
		nTotal = nLengthXY + nLengthIsing;
		numOfThermal = 200 * int(double(lat_.getNumOfSites()) / double(nTotal));
		numOfThermal = std::max(200, numOfThermal);
		double ratio = (nLengthXY / nLengthIsing);
		outFile_ << "Worm Length XY MA - " << nLengthXY
				<< " Worm Length Ising MA - " << nLengthIsing << " Total MA - "
				<< nTotal << " numOfThermal " << numOfThermal << " Ratio "
				<< ratio << " Iter " << numIter << endl;
		numIter++;
	}
	wormLengthXY_ = acc();
	wormLengthIsing_ = acc();
	acc lastTenAv;
	int count = 0, countten = 0;
	for (int i = 0; i < 200; i++) {
		while (ba::count(wormLengthXY_) < 1000
				&& ba::count(wormLengthIsing_) < 1000)
			nextMove();
		if (i > 150)
			lastTenAv(ba::mean(wormLengthXY_) + ba::mean(wormLengthIsing_));
		outFile_ << " Mean - "
				<< ba::mean(wormLengthXY_) + ba::mean(wormLengthIsing_) << endl;
		wormLengthXY_ = acc();
		wormLengthIsing_ = acc();
	}
	int nLength = ba::mean(lastTenAv);
	outFile_ << nLength;
	nTotal = nLength;
	double wormLength = nTotal;
	int closedRate = std::max(10,
			int(
					double(qmcParameters_.rateFactor.getVal())
							* double(lat_.getNumOfSites()) / wormLength));
	measures_.updateClosedRate(closedRate);
	outFile_ << endl << "Closed Rate - " << closedRate << endl;
}

void OnModel::dump(bool final) {
	string outPath = qmcParameters_.outPath.getVal();
	measures_.dump(outPath, final);
	lat_.dump(outPath, final, &outFile_);
	eng_.dump(outPath);
	string fileName = outPath + "/sim_params.xml";
// Write to XML
	pugi::xml_document doc;
	pugi::xml_node rootNode = doc.append_child("SIM_PARAMS");
// Parameters
	pugi::xml_node node = rootNode.append_child("CLOSED_RATE");
	node.append_child(pugi::node_pcdata).set_value(
			(boost::lexical_cast<string>(measures_.getClosedRate())).c_str());
	doc.save_file(fileName.c_str());
	outFile_ << "Finished Dump - " << endl;
//	doc.save(std::cout);
//	cout << endl;
}

void OnModel::read() {
	outFile_ << "Starting Reading From File" << endl;
	string outPath = qmcParameters_.outPath.getVal();
	measures_.read(outPath);
	outFile_ << "Measures Loaded" << endl;
	lat_.read(outPath, true);
	int totalBonds;
	try {
		bool testConf = lat_.testConfiguration(totalBonds);
		outFile_ << "Total Bonds Calculated - " << totalBonds
				<< " total Bonds previous sim - " << lat_.getTotalBonds()
				<< endl;
		if (totalBonds != lat_.getTotalBonds())
			throw "Total Bonds after Read don't match";
		if (!testConf)
			throw "Bad Load of Configuration";
	} catch (string &e) {
		outFile_ << e;
		exit(1);
	}
	outFile_ << "Lattice Loaded" << endl;
	eng_.read(outPath);
	string fileName = outPath + "/sim_params.xml";
// Read From XML
	pugi::xml_document doc;
	try {
		pugi::xml_parse_result result = doc.load_file(fileName.c_str());
		if (!result)
			throw result;
	} catch (pugi::xml_parse_result &result) {
		std::cout << "Error description: " << result.description() << "\n";
		exit(1);
	}
	pugi::xml_node rootNode = doc.child("SIM_PARAMS");
	measures_.updateClosedRate(readXmlData<int>(rootNode, "CLOSED_RATE"));

}
